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Abstract — For nonrigid image registration, matching the par- 
ticular structures (or the outliers) that have missing correspon- 
dence and/or local large deformations, can be more difficult than 
matching the common structures with small deformations in the 
two images. Most existing works depend heavily on the outlier 
segmentation to remove the outlier effect in the registration. 
Moreover, these works do not handle simultaneously the missing 
correspondences and local large deformations. In this paper, 
we defined the nonrigid image registration as a local adaptive 
kernel regression which locally reconstruct the moving image's 
dense deformation vectors from the sparse deformation vectors 
in the multi-resolution block matching. The kernel function 
of the kernel regression adapts its shape and orientation to 
the reference image's structure to gather more deformation 
vector samples of the same structure for the iterative regression 
computation, whereby the moving image's local deformations 
could be compliant with the reference image's local structures. To 
estimate the local deformations around the outliers, we use joint 
saliency map that highlights the corresponding saliency struc- 
tures (called Joint Saliency Structures, JSSs) in the two images 
to guide the dense deformation reconstruction by emphasizing 
those JSSs' sparse deformation vectors in the kernel regression. 
The experimental results demonstrate that by using local JSS 
adaptive kernel regression, the proposed method achieves almost 
the best performance in alignment of all challenging image pairs 
with outlier structures compared with other five state-of-the-art 
nonrigid registration algorithms. 

Index Terms — nonrigid registration, outliers, missing corre- 
spondence, local large deformation, local model, local similarity, 
local structure adaptivity, kernel regression, joint saliency map 
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NOnrigid image registration has attracted increasing at- 
tention at motion tracking, change detection, image 
segmentation and surface reconstruction in the computer vi- 
sion and pattern recognition for the last decade |Q]-|]3]. The 
objective of nonrigid image registration is determining the 
local transformations that align every structures (or features) 
in one (moving) image with the corresponding structures in 
the another (reference) image. However, owing to the image 
content changes over a period of time and the different 
physical mechanisms of multimodal imaging sensors, some 
local structures presented in one image appear partially or even 
disappear completely in another image. Such local structures 
without one-to-one correspondence are closely intertwined 
with the strutures' local large deformations in the nonrigid 
image registration. The missing correspondences and local 
large deformations of structures are callled outliers in this 
paper. At present, obtaining an accurate and robust nonrigid 
image registration is still a challenging task for matching 
these local outlier strucutres with missing correspondences 
and/or the local large deformations. In the computer vision 
community, the nonrigid registration with the local large 
deformations is also knows as the large displacement optical 
flow problem [4|. 

Generally, outliers mean the extreme observations substan- 
tially different from all other ones in the real data. For nonrigid 
image registration, the missing correspondences and local 
large deformations introduce the extreme observations which 
appear as the extreme local geometric and intensity differences 
between the two images to be registed. In geometric morphol- 
ogy, these outliers exhibit some large and complex structural 
discrepancies in the locally varying spatial contexts where the 
underlying different local structures could deform in substan- 
tially different ways. The desired nonrigid image registration 
should be able to match the moving image's local structures 
to the corresponding reference image's structures from these 
various local differences. Therefore, the local regression mod- 
els |5|[6|[35| that are adept at handling the locally varying 
differences are necessary to account for these outliers, and it 
provides the rationale behind this work. Recently, Gerogiannis 
et a/.Q justify that the local Bayesian regression model are 
favorable to model local registration transformation compared 
with other interpolation based registration techniques. 

Over the last several decades, many relevant works were 
proposed to match the local structures (or features) of the 
two images by minimizing the differences between the two 
images, which include feature-based, intensity-based, and 
hybrid methods fDEl , Feature-based ||8]||9] approaches are 
local model based registration methods because they always 
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use the local model of sparse image representation to select 
some corresponding features in the two images, and then 
directly match the local features by finding a geometrical 
transformation from these sparse feature correspondences. 
With the recent development of local feature detectors and 
descriptors (such as SIFT) [4|[10|, the feature-based methods 
have developed into the hybrid methods to use integrated 
local feature signatures [11 1 in characterizing each voxel in 
the two images. Nevertheless, the computation of registration 
transformation is very sensitive to the ambiguity in finding 
local feature correspondences in different local spatial contexts 
with outliers. Moreover, most feature-based methods have not 
solved the outlier issues that have both missing correspon- 
dences and local large deformations. 

By using the complete image data to recover dense cor- 
respondences at pixel-level precision, most intensity-based 
nonrigid registration approaches are regarded as global model 
based methods that are often formulated as global energy 
minimization problems with the energy being composed of an 
regularization term and a similarity term llTZllfPJl . The relative 
weight of similarity term and regularization term can cause the 
well-known trade-off between the registration accuracy and the 
smoothness of the deformation field [14|. In the presence of 
outliers, the accurate and plausible local structure matching 
does not exist using whole-intensity driven transformation 
model. The relative spatial regularization can either cause 
non-smooth and implausible distortion in these outlier regions 
or introduce over-smooth and inaccurate mapping artifacts 
between the whole images. This outlier problem can be solved 
by using a locally varying weight between regularization and 
image similarity [15|-[21|, creating artificial correspondences 
[22|-|27|, using cost-function masking [28|-|30|, or develop- 
ing SIFT flow for large displacement [4 |. These approaches are 
largely dependent on the segmentation of outlier regions and 
give no full consideration to both the missing correspondences 
and local large deformations. 

By successfully handling the locally varying differences 
in pattern recognition and machine learning, local kernel 
regression Il6l ll36l (or nonparametric regression) is regarded 
as an ideal local regression model to effectively reconstruct 
the desired local signal while suppressing the outlier and noise 
effects. The normalized convolution ll36l used by Suarez et al. 
Il37lll38l was the first application of local kernel regression in 
estimating dense deformation fields from sparse deformation 
fields. More recently, two works [39|[40| also utilized ker- 
nel regression to estimate registration transformations. Xing 
and Qiu [42] proposed intensity -based nonparametric local 
smoothing to estimate local mapping transformation in a 
neighborhood. However, these methods do not consider the 
outlier problems in image registration. 

To solve the outlier problem, we proposed the joint saliency 
map (JSM) to group the corresponding saliency structures 
(called Joint Saliency Structures, JSSs) in intensity-based sim- 
ilarity measure computation [31 1. The JSM has been proved to 
greatly improve the accuracy and robustness of rigid (|3"T1[|32"1 
and nonrigid [10|[33|[34| image registration with outliers. We 
further think, by reflecting the local structure correspondence, 
the JSM also could guide the local kernel regression for accu- 



rately estimating registration transformations in the nonrigid 
image registration with outliers. 

In this paper, by integrating the JSM into the kernel regres- 
sion's local adaptive estimation, we propose a new JSS adap- 
tive kernel regression to solve the outlier problem in the non- 
rigid registration. Specifically, with a moving window/kernel 
in kernel regression based nonrigid image registration, the 
output dense deformation vectors are locally estimated based 
on the specific weights of the sparse deformation vectors 
within the moving window. In the presence of outliers, the 
weights of the sparse deformation vectors for the outliers 
should be as small as possible to reduce the outlier effect 
on giving a distorted regression of dense deformations, while 
the JSSs and their underlying sparse deformation vectors 
should be emphasized with their weights being kept as high as 
possible to ensure the precision of the regression computation. 
Furthermore, the kernel function adapts its shape [35|,[44|- 
[47 1 and orientation to the reference image's local structure in 
order to gather more deformation vector samples of the same 
structure in the kernel regression, whereby the regression of 
local deformations can be locally compliant with the underly- 
ing local saliency structures without directing the deformation 
across the edges and corners of local structures. 

To the best of our knowledge, this is the first work which 
propose the local JSS adaptive kernel regression to align the 
local structures by locally estimating the dense deformation 
fields of nonrigid image registration in the presence of outliers. 
An important contribution is that we use the JSM to guide local 
adaptive kernel regression for the accurate matching of local 
structures while suppressing outlier effects on the nonrigid 
image registration. The proposed method also makes the 
regression kernel not only locally adapt to the JSSs in the two 
images but also to the reference image's saliency structures. 
These two adaptivities enhance our local structure matching 
algorithm in achieving the accuracy and robustness of nonrigid 
registration of images with missing correspondences and local 
large deformations. The rest of this paper is organized as 
follows. Our algorithm is elaborated in Section 2 followed 
by experimental results in Section 3. The whole paper is 
concluded in Section 4. 

II. Methods 
A. Coarse-to-fine Block Matching Scheme 

The algorithm proposed by us is built upon a coarse-to-fine 
iterative block matching scheme similar to the one in [37|[38|. 
The coarse-to-fine iterative framework is well known to deal 
with large deformation (or large displacement [4|) in nonrigid 
image registration [12|. Fig. 1 displays the whole coarse-to- 
fine iterative framework, where different levels have their own 
resolutions but the same procedure. At each level, the resulted 
global deformation is composed of initial deformation and 
current deformation. The initial deformation is obtained by 
resampling the global deformation from the previous level 
while the current deformation is the result of the current level 
involving the following two stages. 

In the first stage, a discrete and sparse displacement field is 
created by maximizing local similarity measure for every pixel. 
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Reference Moving 
Image A Image B 

_J L_ 

Block Matching in the First Coarse Level 



Block Matching in the Last Fine Level 

I I 

Resulted Registed Moving 

Deformation Field Image B 

Fig. 1: Flowchart of our algorithm within coarse- to-fine frame- 
work. 



The choice of local similarity measure has been studied in 
our previous work [42|, where we employed cross-correlation 
(CC) instead of mutual information (MI) as the similarity 
measure when we matched two images at lower levels of the 
hierarchy so that the problem of Mi's statistical consistency 
1 43 1 could be solved during the hierarchical subdivision. In 
this study, however, we only employ MI as the local similarity 
measure because we have restricted the block size of the 
lowest level of the hierarchy. Although block matching has 
many advantages in obtaining the deformations of an image, 
implementing only this algorithm is still insufficient to avoid 
the irregularity of transformations such as tearing, folding 
and distorting. Therefore, further reconstruction constraint is 
indispensable to be integrated into the registration procedure. 
To this end, a reconstruction procedure using local kernel 
regression is applied to this discrete and sparse displacement 
fields in the second stage. Details of this stage is described 
in Section 2.2 to Section 2.4. After the above-mentioned two 
stages for each level in this coarse-to-fine framework, a smooth 
and dense deformation field is iteratively achieved as the 
global deformation at the last and finest level. 



B. Kernel Regression Based Local Deformation Reconstruc- 
tion 

Inspired by the successful applications in mordern image 
and video deblurring and super-resolution imaging [35 |,[44|- 
[47], we utilize kernel regression to reconstruct the dense 
deformation fields from the discrete and sparse displace- 
ment fields obtained through block matching. Suppose we 
have sparse and irregularly distributed deformation fields 
{yi)Xi}£=i g iven in the form 

y i =T{x i )+e i , x s: ett, i = l,..-,P (1) 

where y^ is a sparse displacement vector (response variable) at 
position (explanatory variable) Xj, T (■) describes the desired 
dense deformation field in the moving windows (kernel) £1 
with independent and identically distributed zero mean noise 
Ej = £ (Xj). In statistics, the function T (■) is treated as a 
regression of y on x, T(x) = £?{y|x}. In this way, the 
reconstruction of nonrigid deformation fields is from the field 
of the regression techniques. 

Suppose the point of interest x to be reconstructed is near 
Xj, then the regression of dense deformation field X'(xj) can 
be approximated by a local Taylor series expansion 

T( Xi ) a T(x) + {VT(x)} T (x s - x) 

+ -(Xj - x) T {Hessian[T(x)]} T (x j; - x) + ■ ■ ■ 
«A>+/?i'(x i -x) 

+ pf vech{(xj - x)(xj - x) T } + • • • 

(2) 

where vech(-) is a half-vectorization operator of a symmetric 
matrix and {/5o, fli,&2, ■ • • , pV} are N + 1 unknown param- 
eters to be estimated. As in the 2D case, x = [xi,X2] T , we 
can easily estimate the unknown parameters as 

A>=T(x) 

0T(x) dT (x) T 

Pi = — k , — ^ 

ox i dx 2 

1 d 2 T(x) d 2 T (x) ^ d 2 T (x) T W 

2 dx\ ' dx\dxi ' dx\ 

Since we have known the discrete displacement vectors 
{y-t}£=i> we can compute {P n }n=o by finding the optimum 
solution of the following weighted least squares problem 

p 

SR F? ,l][yi-A)-pT(xi-x) 
{ftAA '- } ,=i (4) 

- pfvecri{(xj - x)(Xj - x) T } 

where Kh{-) is a kernel function (see the detail in the next 
section), which not only smoothes the approximation but also 
penalizes distance away from x. 

In addition, if we assume y = [yi,y2,--- ,yp] T , b = 
[/3 ,pT,--- ,I3n] T , and K = diag[^(xi - x),i^(x 2 - 
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x), • • • , Kh(x p — x)], then we can rewrite the optimization 
problem in a matrix form 



b =argmin (y - Xb) J K(y - Xb) 

b 



(5) 



with 



X = 



(I ( Xl -X) 

1 (x 2 - x) 



vech T {(xi — x)(xi — x) T } 
vech T {(x 2 — x)(x 2 — x) T } 



\l (xp — x) vech T {(xp — x)(xp — x) T } 
and the least-squares estimation solution can be expressed as 

b = (X T KX)- 1 X T Ky (6) 

Because the zero-order Taylor series expansion known as 
the Nadaray-Watson estimator is sufficient to reconstruct the 
displacement vectors, the estimation of the deformation field 
at x has the form 



T(x) = p 



(7) 



Ei=i K H (*t - x) 
Since images have outliers, it is reasonable to consider un- 
certainty for each pixel [4-8 1 . Therefore, we add a weight (or 
certainty) function Ci to equation (7) 



T(x) = ft, = 



T,i=i k h(xi ~x) ■ (yj - a) 

Y^=i k h{m - x) • Ci 
K (g> (y • c) 
K <x> c 



(8) 



The last line of equation (8) can also be expressed in the form 
of normalized convolution ||36l . where <g> denotes convolution 
operation. 

Fig. 2. illustrats the discrete displacement vectors recon- 
structed for every pixel in tumor resection region using local 
kernel regression. Because block matching results inherently 
contain incorrect matches, which are exacerbated by the out- 
liers in the tumor resection region. As a result, the conflicts 
between neighboring displacement vectors (see the several red 
circles shown in Fig. 2(a)) widely exist in the discrete displace- 
ment field for the tumor region. Those displacement conflicts 
can easily introduce the topology change of structures, such 
as tearing and distorting. Fortunately, all the displacement 
vector conflicts are removed or smoothed by the local kernel 
regression in Fig. 2(b), where the displacement vectors hav- 
ing opposite directions fully disappear with the displacement 
magnitudes simultaneously being smoothed. Next, to match 
local structures in the presence of outliers, we design structural 
adaptive kernel functions and robust weighing schemes for 
the moving windows/kernels to further boost the accuracy and 
robustness of the local deformation reconstruction. 

C. Local structure-adaptive kernel function 

As a crucial component of local kernel regression, the 
shape of the kernel function (or moving window) determines 
the spatial distribution of samples which are gathered for 
the quality of the locally reconstructed signal. In principle, 
isotropic Gaussian kernels are mostly used as kernel functions 










T - 






J 





Fig. 2: Application of kernel regression in deformation recon- 
struction, (a) Discrete displacement field, (b) Reconstructed defor- 
mation field using kernel regression. 



in nonparametric regression. However, traditional isotropic 
Gaussian kernels are insufficient to cover more samples of 
the same modality along some specific orientations in signal 
reconstruction. Besides, Gaussian kernels' fixed scales and 
orientations can neither detect nor enhance edge structures. 
These factors easily cause diffusion across lines or edges of an 
image. To remedy these restrictions, Pham et al. proposed an 
anisotropic kernel function to adapt its shape to the density of 
sampling |46|. Afterwards, Takeda et al. B71 utilized gradient 
co variance matrices to construct steering kernels, which have 
been proved to possess the ability to capture the edges of an 
image and be extremely robust to noise and perturbations of 
the data. 

In reconstruction of dense deformation field, the local 
kernel function extends along local structure orientation in 
the reference image so that it can gather more samples 
of sparse deformation vectors that correspond to the same 
saliency structures in the reference image. Besides, the kernel 
function contracts in the normal orientations of local saliency 
structures to prevent deformation diffusion across the edges 
between the different structures. Based on these schemes, it 
can effectively reduce the risk of changing the topology of the 
local saliency structures in the kernel regression. Considering 
that local structure tensor (LST) represents the anisotropic 
local saliency structure information of an image[48|[49|, we 
therefore design local structure-adaptive kernel functions using 
the LST information in the reference image. 

Before designing a local structure-adaptive kernel function 
at a certain pixel, we must compute LST in advance to grasp 
the local saliency structure around that pixel in 2D image. We 
assume that 7(x) denotes the intensity value at point x(x, y), 
the gradient information are expressed as I x = |j, I y = 
and V/(x) = [I x I y ] T indicates the local gradient vectors, 
then the gradient structure tensor (GST) in 2D case can be 
described as 

GST(x) = V J(x) V/(x) T = ( j\ ) < 9 ) 



Ixly 



Generally, to integrate the surrounding structural information 
from neighborhoods, the GST is smoothed by a Gaussian filter 
to derive the LST. The scale a of the Gaussian filter is half the 
filter window size, namely a = 1.5, because the size of the 
filter window in our experiments is 3 x 3 pixels (or 3 x 3 x 3 
voxels for the 3D image). Therefore, we derive the LST from 
GST and perform a principal component analysis of the LST 
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as follows: 

LST(x) = G a * GST(x) = A„uu T + A„vv T (10) 

where {A M ,A,;} (X u > X v ) are the eigenvalues with the 
corresponding eigenvectors being {u, v}. These eigenvectors 
contain the information about the local orientation distribution, 
in which pixel values change fastest along u direction while 
slowest along v. Moreover, eigenvalues reflect both intensity 
variation in each direction and morphological information of 
a local region. For example in 2D image, X u w A„ sa 
corresponds to a homogenous region with no measurable 
structure, X u 3> A„ w describes linear structure while 
X u > A„ 3> appears at corner structure. 

With the above-mentioned LST computation in mind, we 
design the kernel functions to meet the following requirement: 
for regions containing obvious structural information such as 
histologic margins in a medical image, the kernels should be 
anisotropic with the regression computation being enhanced 
just along the main orientation and being suppressed along 
other orientations; for homogeneous regions without distinct 
structural information, the kernel should be isotropic with the 
regression computation being equal in all direction. We assume 
that xo denotes a central position in 2D image, x is a position 
in its neighborhood, {u, v} are the eigenvectors of LST(kq) 
and {A„,A,;} are the corresponding eigenvalues of {u, v}, 
thus a local structure-adaptive Gaussian kernel in 2D case is 
designed as 

a (x,x ) = — — exp 

ZTT<7 u a v 

d u = (d, u), d„ = (d, v 

where d is the vector from Xo to x, {d u ,d„} are the 
projections of the vector d on {u, v}, and the directional 
scales of the Gaussian kernel {a u , a v } are determined by the 
anisotropy A as follows 

a a+A n ^ 
c-u = — ; — rCc, a v = a c (12) 

a + A a 

where A = (X u — X V )/(X U + X v ) . The two directional scales 
of the Gaussian kernel can be adjusted by the parameter 
a > and the local scale a c . Specifically, a determines 
the eccentricity of the kernel while a c affects the number 
of discrete vectors that contribute to the reconstruction of 
continuous deformation vectors. For the sake of simplicity, 
the local scale is set to half the neighborhood window size for 
each kernel according to our experience. We also set a = 0.5 
and <j c = 1.5 because we utilize a 3 x 3 pixel neighborhood 
window in our experiments. Similarily, the methodolgy of 3D 
LST and anisotropic kernel function computations is presented 
as an Appendix to this paper. 

Two kernel functions for two distinct image structures are 
displayed at the doll images in the Fig. 3, where the crosses 
indicate two different kinds of typical image structures (blue 
cross for border and red cross for homogeneous region). In 
Fig. 3, the two pairs of orthogonal vectors indicate the main 
axes of their corresponding Gaussian kernels (blue cross at 
Fig. 3(c) and red cross at Fig. 3(d)). The length of the vector 
is determined by the scale in the direction of the vector. 









2a u 


2<tJ_ 


(11) 


= X 


- x 






+ 












Fig. 3: Gaussian kernels designed for different image local 
structures, (a) Two labeled positions (red cross and blue one), (b) 
The scales and orientations of Gaussian Kernels in corresponding 
positions, (c) Gaussian kernel for the region with blue cross, (d) 
Gaussian kernel for the region with red cross. 
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Fig. 4: Comparison between using isotropic kernels and using 
local structure-adaptive kernels in kernel regression based defor- 
mation reconstruction, (a)-(b) The reference and moving images, 
(c)-(d) Registered images using isotopic kernels and using local 
structure-adaptive kernels in kernel regression, (e)-(f) Local displace- 
ment vector fields for (c)-(d), (g)-(h) Global deformation mesh fields 
for (c)-(d). 



Fig. 4 illustratively explain why we prefer local structure- 
adaptive kernel to isotropic kernel in our structure-adaptive 
kernel regression. The regions pointed by red arrows are 
small scale structures which have local large deformations. 
The directions of the displacement vectors (spaced every 5 
pixels) in these small structures are inevitably conflicted with 
those of the large deformations of the neighboring structures. 
These defomration conflicts that are introduced by the opposite 
displacement vectors can easily cause tearing, folding or 
distorting of the local small scale structures. For example, the 
eyes in Fig. 4(c) display distortion owing to the conflict of the 
deformation directions displayed in Fig. 4(e). Comparatively, 
our local structure-adaptive kernels suppress the contributions 
of the displacement vectors which are not consistent in the 
structure orientation, the deformation conflicts are therefore 
removed in Fig. 4(f) with no distortion existing in the eyes at 
Fig. 4(d). 

Fig. 4(h) demonstrates the overview of mesh deformation 
(10 pixels vertex spacing) in the local structure-adaptive kernel 
regression, which can produce the smooth adaptivity of local 
meshes deformation to the local anisotropic structures. Ob- 
viously, this local structure-adaptive kernel regression obtain 
smooth mesh boundaries which are consistent with the local 
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structures' boundaries. However, isotropic kernel in kernel 
regression easily produce irregularly deformed local meshes 
which are not adaptive to the local structures, so that it is 
very difficult to identify boundaries of local structures from 
the non-smooth mesh deformation in Fig. 4(g). 

D. Robust Weight Mechanism using JSM 

In original kernel regression, the weight c(x), between 
and 1, specifies the reliability (or certainty) at x for the 
local estimation in a moving window/kernel. In modern local 
regression model, the weight function always describes the 
spatial dependence in the locally varying special context. 
With locally data-dependent weights, recently popularized and 
very effective image filter are developed in image and video 
processing field (6). For example, the Gaussian function of a 
residual error with an acceptable range of error a r 11461 severs 
as a new weight function in the neighborhood of x 

|J(x)- J(x,xp)|% 
c x, x ) = cxp( — (13) 

where 7(x) denotes a measured intensity at position x, and 
7(x, xo) is an estimated intensity at x using an initial poly- 
nomial expansion at Xo. 

Since the regions with salient structure information have 
real influence over locallly adaptive image processing based 
on kernel regression, E. Suarez et al. [37] proposed a weight 
scheme in the kernel regression of registration transformations 
by using the scalar measure of a local structure in reference 
image. However, for nonrigid image registration, the salient 
structural regions in reference image may introduce non- 
corresponding salient structural regions at same locations in 
moving image. Therefore, the method proposed by E. Suarez 
et al. l37l is most likely to assign wrong high weights to the 
salient structures that have missing correspondences. 

To avoid the above-mentioned mis-assignment and mini- 
mize the ourlier impact on the reconstruction of deformation 
field, we propose a robust weight mechanism by simul- 
taneously considering the matching degree of local salient 
structures in the overlapping parts of the two images. In our 
previous work J3T), it has been proved that the application 
of JSM is effective in tackling registration problems with 
outliers by emphatically grouping the JSSs into intensity-based 
similarity measure. Continuing the success of JSM, we further 
deploy the concept of JSM into our robust weight mechanism 
and pay more attention to the JSSs in the two images. Those 
JSSs and their incurring deformation should be emphatically 
treated in the local kernel regression to reconstruct the dense 
deformation fields from the sparse deformation fields that 
contain outliers. 

With the general JSM-based weight scheme in mind, we 
should first construct a saliency map to indicate the local 
salient structure distribution in each image. Since LSTs men- 
tioned in Section 2.3 contain sufficient structural informa- 
tion in a region, we can reasonably use them to reflect the 
saliency map of an image. Inspired by the center-surround 
mechanism [50] which has defined the intensity-contrast-based 
visual saliency map, we define a saliency operator based 



on the contrast among neighboring structure tensors. This 
contrast emphasizes the dissimilarity or discrepancy between 
neighboring local structure tensors. Specifically, for a given 
point xo and its neighborhood ft, the saliency value 5(xo) at 
xo in a salient map can be computed through 

5(x ) = avg^ \\LST(x) - LST(xo)\\ D (14) 

xGfi 

where || • defines a distance metric describing the dissim- 
ilarity between two LSTs, which is detailed in the follow- 
ing section. The operator avg computes the average of the 
dissimilarities within the neighborhood SI of Xq. Traditional 
tensor similarity measures such as fractional anisotropy (FA) 
and cosine similarity measure are not appropriate for defining 
tensor-based saliency operator because they only compare 
either scales or orientations of two tensors. Fortunately, some 
improved tensor similarity measure computing both scale 
information and orientation information have been reported. 
In BTl . H. Zhang et al. introduced diffusion tensor metric, 
which is defined as 



||Ti - T 2 \\ L = ^(\\Ti - T 2 \\l + -Tr*{T x - T a )) (15) 

where \\T X - T 2 \\c = V Tr ( T i - T z) 2 is the Euclidean 
distance between two tensors {T\, T 2 }, Tr is the operator for 
computing the trace of matrix. Afterwards, H. Zhang et al. 
[52 1 modified equation (15) to only focus on the anisotropic 
components between two tensors. The modified equation that 
is adopted at equation (14) can be expressed as 



\\Ti T 2 \\ D = \J^{\\Ti - T 2 \\l - -Tr 2 (Ti - T 2 )) (16) 

In a saliency map, the saliency values are general rep- 
resentation of the local structure distribution in an image. 
Low saliency values always appear in the homogeneous and 
background regions while high saliency values are assigned to 
the edge regions of saliency structures owing to the highligted 
contrast among neighboring LSTs in these regions. After 
the two normalized salient maps were achieved to indicate 
the local structure distribution, JSM is ready to describe the 
matching degree between the two saliency maps at every pixel 
pair in the overlapping regions of the two images. Given a 
point Xr in the reference image and its corresponding point 
Xjvf in the moving image after initial transformation, their 
joint-saliency value in a JSM is defined as 

J5(x fl ,x M ) 

A ■ B 

= min{5 , fi (x fl ),5 , M (x M )} ,, , — r rr— 

B + \\LST(x R ) - LST{tc m )\\d 

(17) 

where {Sr(-), Sm(-)} denote the saliency values in the 
saliency maps of the reference and the moving images. The 
empirical parameter A and B are used to normalize the JSM 
values into a final value between and 1 . In our experiments, 
A = 10 and B = § max(||iS , T(x i? ) - LST{x M )\\ D . It 
should be noted that it may introduce a situation that both 
of the corresponding pixels are assigned high saliency values 
in the saliency maps, while their local variations of gradient 
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Fig. 5: JSM Examples with colour scale representing different 
joint saliency values. Each column shows one case. The reference 
images and the moving ones are shown at the top and middle rows. 
Their JSMs are displayed at the bottom rows where the red regions 
correspond to the higher joint saliency values while the blue regions 
correspond to the lower joint saliency values. 
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Fig. 6: The reference and the moving images and their gradient 
magnitude, largest eigenvalue profiles for GST and LST, and 
the JSM magnitude, (a)-(b) The reference and moving images, (c)- 
(e) Gradient magnitude profiles, largest eigenvalue profiles of GST 
and largest eigenvalue profiles of LST of the red line in (a), (f)-(h) 
Gradient magnitude profiles, largest eigenvalue profiles of GST and 
largest eigenvalue profiles of LST of the red line in (b), (i)-(j) Saliency 
value profiles of the red lines in (a)-(b), (k) JSM value profiles of 
the red lines in (a)-(b). 



orientations are in fact totally different. To avoid this situation, 
we also consider the dissimilarity measure between LST(xr) 
and LST(xm) at the denominator in equation (17). 

Fig. 5 shows some examples of normalized JSM with the 
colour scale representing different joint saliency values. The 
high joint saliency values being represented by red colour 
suggest that the underlying pixel pairs come from the JSSs. 
On the contrary, the regions with low JSM values are rendered 
in blue colour, which indicates that the underlying pixel pair 
originates from either homogeneous regions or outlier regions. 
The discrete displacement vectors in these red JSS regions 
are expected to contribute more to the kernel regression than 
the blue regions having low JSM values, this scheme is 
therefore called JSS adaptive kernel regression for nonrigid 
image registration. 

The JSM in our study mainly responds to the corresponding 
high-gradient edge pixels. However, it does not simply high- 




Fig. 7: Comparison between using JSM-based robust weight 
mechanism and not using it in our method, (a)-(b) The reference 
and moving image, (c)-(d) Registered images without and with using 
a JSM-based robust weight mechanism, (e)-(f) Local displacement 
vector fields of (c)-(d), (g)-(h) Global deformation mesh of (c)-(d). 



light the common image gradients in the two images. Fig. 6 
presents the image gradient magnitude, the largest eigenvalues 
of GSTs and LSTs, the saliency values and the JSM values 
profiles of the same red line at the two images (Fig. 6(a)- 
(b)). For easy comparison, the range of ordinates in Fig. 
6(c)-(k) are bound to [0,1]. As shown in Fig. 6, the image 
gradient features in Fig. 6(c) and Fig 6(f) are very sensitive 
to noise and do not agree with each other at each overlapping 
leocation. The noise sensitivity are gradually reduced by using 
the GSTs (Fig. 6(d) and Fig. 6(g)) and the LSTs (Fig. 6(e) 
and Fig. 6(h)). The saliency values of the two images in 
Fig. 6(i)-(j) are robust to noise due to their computing the 
regional contrast of LSTs through equation (14). Moreover, 
the structural image information in a large region is also 
comprehensively considered according to equation (14). As 
a result, the JSM values (Fig. 6(k)) computed through the 
saliency values can accurately preserve the JSSs in larger 
capture range with smaller variability than the image gradients. 

Because of the outliers introduced by missing correspon- 
dence, local large deformation and incorrect block matching, 
the dense deformation fields can not be interpolated form 
the sparse displacement vectors in block matching. The JSS 
adaptive kernel regression is then used to reconstruct the 
dense deformation fields from the sparse displacement vec- 
tors, i.e., smooth the impact of outliers on the deformation 
reconstruction. Due to the expected deformation in the outlier 
region being consistent with its neighboring deformations, the 
JSM in the neighboring regions is used to assign different 
weights to the different displacements of the neighboring 
structures, only those neighboring deformations with high JSM 
values indicating the consistency in structure orientations are 
given high weights in kernel regression based deformation 
reconstruction. 

Fig. 7 illustrates the improvements on the deformation 
field reconstruction after introducing the JSM-based robust 
weight mechanism in the local JSS adaptive kernel regres- 
sion. The regions pointed by red arrows in Fig. 7(c)-(d) are 
outlier regions with missing correspondences and local large 
deformations. Without JSM-based robust weight mechanism, 
the converged displacement vectors (5 pixels spacing) from 
conflicting directions (See Fig. 7(e)) in the outlier region 
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spread the distortion effect into the eye region (see Fig. 7(c)). 
Due to the JSM-based robust weight mechanism introducing 
weighted smoothing effect on the magnitude and direction 
of displacement vectors (see Fig. 7(f)), the moving image's 
eye distortion is removed with the outlier region's structure 
deformations being simultaneously matched to those of the 
reference image (see Fig. 7(f)). Compared with the defor- 
mation meshes (10 pixels vertex spacing) in Fig. 7(g), the 
deformation meshes at Fig. 7(h) display the overall smoothness 
improvement for the structural deformations at the outlier 
regions. 

III. Experimental Results 

In comparison with several state-of-the-art nonrigid reg- 
istration approaches to validate the proposed algorithm on 
some challenging images with missing correspondences and 
local large deformations, we choose Advanced Normalized 
Tools (ANTsf] with elastic transformation and MI (AMI) (53] 
due to the ANTs being placed the first in the EMPIRE 10 
challenge [56 1 evaluating 34 total registration algorithms. We 
also include Diffeomorphic Demons with Diffusion-like reg- 
ularization (DDDjl 154), fast B-Spline with MI (BMlfl AMI 
with cost-function Masking (AMM) and Large displacement 
Optical Flow (LOF) algorithms^ [4]. The parameters of our 
methods are: the number of pyramid levels is 5; the local sim- 
ilarity measure is mutual information. We set the parameters 
of AMI and AMM as: the histogram bin size is 32; the number 
of pyramid levels is 3; the iterations are set to 100 x 100 x 10; 
the gradient step is 10; the default regularization is Gaussian 
filtering with a sigma of 3. The parameters of DDD method 
are set as follows: the variance of smoothing kernels is 2; the 
step scale is 1; the number of pyramid levels is 5; the number 
of iterations is 100. The parameters of the BMI method are 
selected as: the number of iterations is set to 100; the grid 
size is 15; the histogram bin size is 32; the spatial sample is 
50000; the maximum deformation is 20. We choose the default 
parameters of the LOF method as: the tuning parameters for 
regularity constraint, the point correspondence constraint and 
the contraint on the gradient are 30, 300 and 5; downsampling 
factor is 0.95; the numbers of outer iterations and inner 
iterations are set to 10 and 5. With those parameters all the 
methods mentioned above achieve their best performances. 

To evaluate the performance of the six competing methods, 
both registration accuracy and efficiency are estimated during 
the assessment. Validating the registration accuracy of binary 
image pairs is easy to recognize the badly-aligned regions 
by using the difference image between the reference and 
the registed moving images. However, the evaluation based 
on the difference images is not reliable for grayscale image 
registration (551 . Due to the registration errors measured at 
densely distributed landmarks being considered as the standard 
for evaluating registration accuracy of grayscales images [55 1, 
we manually selected a large number of appropriate landmark 

'http://www.picsl.upenn.edu/ANTs 

2 http ://mipav. cit. nih.gov 

3 http ://w ww. sheer, org 

4 http://www.seas.upenn.edurkatef/LDOF.html 
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Fig. 8: Chinese character image registration, (a)-(b) The reference 
and moving images (c) the proposed method, (d) AMI, (e) DDD, (f) 
BMI, (g) AMM, (h) LOF. 
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Fig. 9: Difference images from the six Chinese character reg- 
istration results, (a) the proposed method, (b) AMI, (c)DDD, (d) 
BMI, (e) AMM, (f) LOF. 



pairs in the reference and the registered moving images to 
accurately evaluate the registration accuracy of the grayscale 
images. The selected landmarks fully exclude outlier features 
but identify some corresponding small scale saliency structures 
(having local large deformations) around the outlier regions. 
Therefore, the matching accuracies of six registration methods 
for the local structures with outliers are fully assessed by using 
the Mean Registration Error (MRE) and Standard Deviation 
(SD) in pixels between these selected landmarks. 

The reference and the moving binary images in the first 
experiment are two similar Chinese characters at Fig. 8(a)- 
(b). The two strokes in the left part of the reference image 
correspond to the upper left and lower left ones in the 
moving images, while the middle stroke (outlier) in the left 
part of moving image has missing correspondence (see the 
red rectangles in the Fig. 8). Moveover, the triangular and 
the rectangular openings at the right part of the character 
are narrowed. The local large deformations are apparent at 
the lower left comer of the rectangular region and the low 
right corner of the triangular region. Besides, the lower left 
stroke from southwest to northeast is lengthened with counter- 
clockwise rotation. Fig. 8(c)-(g) are the registered results of 
our approach, AMI, DDD, BMI, AMM and LOF. The strokes 
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Fig. 10: Mickey image registration, (a)-(b) The reference and 
moving images, (c) the proposed method, (d) AMI, (e) DDD, (f) 
BMI, (g) AMM, (h) LOR 




Fig. 11: Brain tumor resection image registration, (a)-(b) The 
reference and moving images, (c) the proposed method, (d) AMI, (e) 
DDD, (f) BMI, (g) AMM, (h) LOF. 



at the left part and the corner regions of the two openings are 
locally deformed more or less in the six registered images. 
Though the constrained cost-function masking has masked the 
middle stroke in moving image, AMM method perform poorly 
in matching local structures with local large deformations in 
Fig. 8(g). The LOF method changed the topology of moving 
image by fully removing the middle stroke in the left part 
and introducing some artifacts in the character at Fig. 8(h). 
Overall, the proposed method (see Fig. 8(c) has produced the 
best structure deformation among the six registered results. 

The performance of our proposed approach could be clearly 
validated from the difference images (see Fig. 9) of the six 
registration results, where the white regions represent the 
discrepancies between the reference image and the registered 
moving image. Althogh the difference image of LOF method 
has smallest white regions among the six methods, the LOF 
method also introduce distinct artifacts and topology change 
in the registered Chinese character. Without robust scheme 
tackling missing correspondence and local large deformation 
simultaneously, the AMM method (Fig. 9(e)) even performs 
a little worse than the AMI method (Fig. 9(b)). The white 
regions in the difference image (Fig. 9(a)) of our approach 
are the least among the other five results which preserve the 
topology of the middle stroke. It validates that our proposed 
method match well the moving image's local structures into 
the corresponding structures at the reference image. 

The second experiment involves aligning two grayscale 
Mickey images, which have the outlier doctoral hat shown in 
the moving image (Fig. 10(b)). Moreover, the local structure's 
large deformations occur in Mickey's left thumb, left hand, 
right thumb (see the red rectangles in the Fig. 10), right shoe 
and right button in Mickey's belly. Consequently, validating 
the performance of registration methods is largely dependent 
on the deformed results of these structures. Fig. 10(c)-(h) 
show the registered results of the six methods, where our 
proposed approach outperforms the other five methods by 
perfectly deforming those local structures to desired positions. 
In contrast, the morphologies of Mickey's left hand in Fig. 
10(d) and (g) are changed by AMI and AMM methods. The 
Mickey's left thumb, right shoe and right button in Fig. 10(e) 
have not changed by DDD method. The Mickey's left thumb 
has not deformed and the left palm become thicker in Fig. 




Fig. 12: Flower image registration, (a)-(b) The reference and 
moving images, (c) the proposed method, (d) AMI, (e) DDD, (f) 
BMI, (g) AMM, (h) LOF. 



10(f) after BMI registration. In this case, though with the 
constrained cost-function masking for the doctoral hat, AMM 
method has no improvement in matching local structures (Fig. 
10(g)). The LOF method introduces an undesired artifact in the 
right thumb of the Mickey (see the left red rectangle in the 
Fig. 10(h)). 

Another grayscale image registration is matching pre- and 
post-operative brain tumor images. Brain tissue severely sup- 
pressed by tumor in the preoperative image (Fig. 11(b)) 
expands after tumor resection. Tumor resection not only 
brings missing correspondences into the tumor region of the 
post-operative image (Fig. 11(a)) but also incurs local large 
deformations that are caused by brain shift. A successful 
registered result of this case should properly deform pre- 
operative brain tissue according to the post-operative image 
regardless of tumor resection. Visual inspection has revealed 
that the proposed, AMI, AMM and LOF methods (see Fig. 
ll(c)-(d),(g)) apparently perform better than the DDD and 
BMI methods (see Fig. ll(e)-(f)) because the local brain 
deformation resulted from the latter two methods is either 
insufficient or somewhat excessive. This visual valuation is 
further confirmed by validating landmark-based registration 
error in the following section. 

A more challenging flower image registration is displayed 
in Fig. 12, where the small scale stamen filament in the right 
part of the image is largely deformed while some buds behind 
the stamen filament of the moving image being disappeared 
in the reference image (see the top red rectangles in the 
Fig. 12). Matching small scale structures that have large 
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TABLE I: Registration errors (Mean+SD) of the six methods for 
the three grayscale image registrations 

proposed AMI ODD — BMl AMM EOF 

1.27±3.09 1.56±3.38 6.07±5.93 3.23±5.51 — 12.03±14.80 3.55±3.21 
0.97±1.91 0.95±1.48 1.60±3.08 1.15±1.89 1.16±1.99 1.97±1.93 
1.14±2.96 1.69±3.49 8.38±8.82 4.42±5.65 1.55±3.49 2.72±2.35 



deformations and missing correspondence is very difficult for 
nonrigid image registration. In the case of Fig. 12, an ideal 
registration algorithm should accurately match not only the 
large petal (see the bottom red rectangles) but also the small 
scale stamen filament while simultaneously keeping reasonable 
local deformation consistency in the spatial context. Among 
these six methods, only the proposed approach aligns not only 
the small scale stamen filament but also the large scale petal 
by computing their resonable and accurate deformations. The 
AMI (Fig. 12(d)) and AMM (Fig. 12(g)) methods achieve 
better registration performances than the DDD and BMI meth- 
ods (Fig. 12(e)-(f)). The LOF (Fig. 12(h)) method gets good 
deformations for most structures, but also apparently intruduce 
some artifacts in the petals of the registered moving image. 



The MRE and SD of the manually selected landmarks for 
the six methods in the three grayscale image registration are 
listed in Tab. 1. The proposed method for the Mickey image 
achieves the smallest registration error of 1.27 ± 3.09 pixels 
while the registration errors of AMI, DDD, BMI, AMM and 
LOF methods are greater than or equal to 1.56 ± 3.38 pixels. 
Compared with other method, the proposed method and AMI 
have achieved sub-pixel registration accuracy for the brain tu- 
mor resection images with the registration errors of 0.97±1.91 
and 0.95 ± 1.48 pixels, respectively. As for the flower image, 
the proposed method gets the smallest registration error of 
1.14 ± 2.96 pixels, while the registration errors of other five 
methods are greater than or equal to 1.69 ± 3.49 pixels. 
In average, the proposed method maintains almost the best 
performance in comparison with other five methods. Although 
the orignial AMI method in the brain tumor resection image 
registration has a slight advantage over the proposed method, 
using cost-function masking makes AMM method worse than 
the proposed method in matching locally deformed structures. 
Simply setting the brain-tumor-region's value to zero by cost- 
function masking is not enough to accurately match local 
salient structures with missing correspondences and local large 
deformations. Due to the LOF method only considering the 
effect of large displacments, its performance is not desired for 
the nonrigid registration with both missing correpondences and 
local large deformations. 

Table 2 lists the computation time needed for the six 
algorithms at the four image registration, where the image 
resolution of case 1-3 is 372x392 pixels, and that of case 4 
is 384x288 pixels. All the six methods are operated on a PC 
of Pentium(R) Dual-Core 3.20 GHz CPU with 2 GB memory. 



TABLE II: Computation runtime in seconds for the six registration 
methods (Pentium(R) Dual-Core 3.20 GHz RAM 2.0 GB) 



Cases 


proposed 


AMI 


DDD 


BMI 


AMM 


LOF 


1 


48 


38 


12 


31 


38 


112 


2 


49 


56 


13 


42 


45 


144 


3 


48 


21 


12 


43 


20 


113 


4 


36 


21 


10 


75 


7 


168 



IV. Conclusion 

In this paper, considering the outlier structures that have 
missing correspondences and local large deformations in non- 
rigid image registration, we use JSS adaptive kernel regression 
to reconstruct dense deformation field (for the moving image) 
from the sparse deformation field hierarchically computed 
from multi-resolution block matching. Specifically, the pro- 
posed local JSS adaptive kernel regression implements two 
local adaptivities into the underlying saliency structures in the 
reference images and the JSSs in the two images to accurately 
and robustly match corresponding local structures with missing 
correspondence and local large deformation. 

Nevertheless, there are still some issues that need to be 
further addressed in our future work. First, the local scale 
mentioned in Section 2.3 is set as a constant value for the sake 
of simplicity at each anisotropic Gaussian kernel. As shown 
in some registered results, however, the deformations of some 
small scale structures (such as the eyes in the doll image at Fig. 
8) are affected somewhat badly by the deformations of the sur- 
rounding large scale structures. This situation is caused by the 
unchangeable local scales. Indeed, the local scale (the width of 
kernel) controls the number of discrete displacement vectors 
contributing to the reconstruction of dense deformation vectors 
from sparse displacement vectors. Therefore, the choice of 
local scale significantly affects the final registration result. 
A self-adjustable local scale according to the local structure 
properties is expected to automatically adjust the number of 
discrete displacement vectors participating in the deformation 
reconstruction. To our knowledge, almost no attention has 
been paid on the self-adjustable local scale for nonrigid image 
registration during the last decade. 

Second, we could also improve the proposed algorithm 
by replacing the block matching with other feature-based 
or hybrid nonrigid registration algorithms so that we can 
accurately initialize deformation field estimation and apply 
the proposed method to match multimodal images. As for the 
diffeomorphism of nonrigid image registration [53 1 [54] [56 1, 
the missing correspondences and local large deformations 
make it unrealistic for nonrigid registration method enforcing 
the diffeomorphic local structure matching. However, we could 
initialize our method with some diffeomorphic registration 
algorithms to find the local structures' optimal diffeomorphic 
matching for the subsequent JSS adaptive kernel regression. 
At last, further researches are required to reduce the compu- 
tational cost of our approach. At present, more than half of 
the total running time for our proposed method is spent on 
the local adaptive Gaussian kernel reconstruction and subse- 
quent adaptive kernel regression at every pixel. We should 
design fast method |f57l to estimate discrete local structure- 
adaptive Gaussian kernel and implement structural adaptive 
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kernel regression at every pixel. All these improvements could 
prompt us to make further contributions to the nonrigid image 
registration. 

Appendix A 

3D LOCAL STRUCTURE TENSOR AND ANISOTROPIC 
KERNEL FUNCTION 

The 3D LST can be computed for every point x(x, y, z) in 
3D image as follows. With the gradient being computed as 

7 * = §■ h = §■ 7 * = I and the v/ ( x ) = % h Q T 

indicating the local gradient vectors, the 3D GST is estimated 



GST(x) = V/(x)W(x) 



Ixly 




(18) 



while the smoothed 3D LST being 

LST(x) = G a * GST(x) = A„uu T + X v vv T + A ffi ww T 

(19) 

where {X u , X v , X w } (X u > At, > X w ) are the eigenvalues and 
their corresponding eigenvectors are denoted as {u, v, w}. The 
voxel values change fastest along u direction while slowest 
along w direction. 

To derive the 3D kernel function at xo, we express the 
anisotropic 3D Gaussian function a(x,xo) using the eigen- 
values and corresponding eigenvectors of 3D LST. 



exp 
d, 




(d,u), d„ = (d,v) 



2 d 2 


d 2 

w 


t u la v 


2a w 


w) , d = x 


- x 



(20) 

where d is the vector from Xo to x, {d^d^d^,} are the 
projections of the vector d on {u, v, w}. In order to estimate 
the directional scales of the 3D anisotropic Gaussian kernel 
[58 1, we first compute the anisotropy {a vw ,a uw } of 3D LST 
with 

^ Xv X w ^ X u X w (21) 

An + Xv "+" X w X u — Xv "+" X w 

We further estimate the spatial dependent corner strength 

C(x ) as 

C(x ) - {l-a vw ~a uw )\VI(x )\ 2 (22) 

where |V/(xo)| 2 is the local gradient strength at the point 
xo- Therefore, we obtaint the three directional scales of the 
3D Gaussian kernel as 

_ 0"c(l 0>vw Q"uvj) 



1 + C(xo) 
cr c (l — 2a vw ) 
1 + C(x ) ' 



(23) 



1 + C(x ) 
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